Préparation du jeu de données

On loade la matrice de reads par échantillon et par gènes obtenue par htseq-count. On supprime les dernières lignes qui correspondent aux reads non identifiés comme des gènes ou comme ambigus.

data = read.table("quantif_all_unstranded.txt", h = F)
colnames(data) <- c("Gene", "WT.K1", "WT.K2", "WT.N1", "WT.N2", "WT.K3", "WT.N3")
rownames(data) = data$Gene
data = dplyr::select(data, -c("Gene"))
genes = which(!(grepl("__", rownames(data))))
data = data[genes, ]
data = data[, sort(colnames(data))]
kable(head(data))

Normalisation

In the biological point of view, a gene must be expressed at some minimal level before it is likelyto be translated into a protein or to be biologically important. In addition, the pronounceddiscreteness of these counts interferes with some of the statistical approximations that areused later in the pipeline. These genes should be filtered out prior to further analysis.As a rule of thumb, genes are dropped if they can’t possibly be expressed in all the samplesfor any of the conditions. On met aussi bien à jour la nouvelle taille de librairie.

y <- DGEList(counts = data, group = c(1, 1, 1, 2, 2, 2))
y$samples

      group lib.size norm.factors
WT.K1     1 12527631            1
WT.K2     1  9896571            1
WT.K3     1  3782109            1
WT.N1     2  9246242            1
WT.N2     2  9428460            1
WT.N3     2 11087062            1

keep <- filterByExpr(y)
y <- y[keep, , keep.lib.sizes = FALSE]

RNA-seq provides a measure of the relative abundance of each gene in each RNAsample, but does not provide any measure of the total RNA output on a per-cell basis. This commonly becomes important when a small number of genes are very highly expressed in one sample, but not in another. The highly expressed genes can consume a substantial proportionof the total library size, causing the remaining genes to be under-sampled in that sample.

The calcNormFactorsfunction normalizes for RNA composition by finding a set of scalingfactors for the library sizes that minimize the log-fold changes between the samples for mostgenes. The default method for computing these scale factors uses a trimmed mean of M-values (TMM) between each pair of samples [30]. We call the product of the original library size and the scaling factor the effective library size. The effective library size replaces theoriginal library size in all downsteam analyses. TMM is recommended for most RNA-Seq data where the majority (more than half) of thegenes are believed not differentially expressed between any pair of the sample

y <- calcNormFactors(y)
y$samples

      group lib.size norm.factors
WT.K1     1 12506893    1.0476659
WT.K2     1  9880693    0.9668647
WT.K3     1  3776180    0.9149596
WT.N1     2  9229304    1.0529387
WT.N2     2  9411525    1.0226377
WT.N3     2 11068144    1.0020391

norm <- cpm(y, normalized.lib.sizes = TRUE)
not_norm <- cpm(y, normalized.lib.sizes = FALSE)
rd_genes = sample(rownames(norm), 500)
heatmap(not_norm[rd_genes, ], main = "Sans normalisation")

heatmap(norm[rd_genes, ], main = "Avec normalisation")

The normalization factors of all the libraries multiply to unity. A normalization factor belowone indicates that a small number of high count genes are monopolizing the sequencing,causing the counts for other genes to be lower than would be usual given the library size. Asa result, the library size will be scaled down, analogous to scaling the counts upwards in thatlibrary. Conversely, a factor above one scales up the library size, analogous to downscalingthe counts.

The GC-content of each gene does not change from sample to sample, so it can be expectedto have little effect on differential expression analyses to a first approximation. Recent pub-lications, however, have demonstrated that sample-specific effects for GC-content can bedetected [28, 13]. The EDASeq [28] and cqn [13] packages estimate correction factors thatadjust for sample-specific GC-content effects in a way that is compatible with edgeR.

Like GC-content, gene length does not change from sample to sample, so it can be expectedto have little effect on differential expression analyses. Nevertheless, sample-specific effectsfor gene length have been detected [13], although the evidence is not as strong as for GC-content.

Estimation de la variance et Expression différentielle

Biological Variance is the coefficient of variation with which the (unknown) true abundanceof the gene varies between replicate RNA samples.

The dispersion estimates can be viewed in a BCV plot:

y <- estimateDisp(y)

Design matrix not provided. Switch to the classic mode.

et <- exactTest(y)
best = topTags(et, n = 2000)
best = best[best$table$FDR < 0.05, ]
genes <- rownames(cpm(y, normalized.lib.sizes = TRUE))
DEgenes = rownames(best$table)
plotBCV(y)

plotMD(et)
abline(h = c(-1, 1), col = "blue")
abline(h = c(0), col = "red")

summary(decideTests(et))

         2-1
Down     748
NotSig 14091
Up       928

On a 102 gènes différentiellement exprimés avec un seuil pour la pvalue ajustée de 5%.

The functionplotMDS draws a multi-dimensional scaling plot of the RNA samples in whichdistances correspond to leading log-fold-changes between each pair of RNA samples. Theleading log-fold-change is the average (root-mean-square) of the largest absolute log-fold-changes between each pair of sample

plotMDS(y, top = 100, col = rep(1:2, each = 3))

logcpm <- cpm(y, log = TRUE)
heatmap(logcpm[DEgenes, ])

On peut maintenant chercher à voir les ontologies des gènes sélectionnés :

library(knitr)
mart = useMart(biomart = "plants_mart", host = "plants.ensembl.org", dataset = "athaliana_eg_gene")
results <- getBM(filters = "ensembl_transcript_id", attributes = c("ensembl_transcript_id", 
    "description", "external_gene_name"), values = DEgenes, mart = mart)
rownames(results) = results$ensembl_transcript_id
r = results[DEgenes, ]
r$pval = best$table$PValue
kable(r, caption = "Genes DE et leur GO")